/*
 * topological-hex, a program to compute combinatorial hexahedral meshes.
 *
 * Copyright (C) <2018> <Université catholique de Louvain (UCL), Belgique>
 *
 * List of the contributors to the development, description and complete
 * License: see LICENSE file.
 *
 * This program (topological-hex) is free software:
 * you can redistribute it and/or modify it under the terms
 * of the GNU General Public License as published by the Free
 * Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program (see COPYING file).  If not,
 * see <http://www.gnu.org/licenses/>.
 */

#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <ctype.h>

#include "mesh.h"
#include "errors.h"

static error_code read_vertices(mesh *mesh, FILE *file);
static error_code read_quads(mesh *mesh, FILE *file);
static error_code read_hexes(mesh *mesh, FILE *file);


static error_code reorient_quads(vertex_index *quads,
                                 quad_index num_quads,
                                 vertex_index num_vertices);

void mesh_init(mesh *mesh) {
  mesh->num_vertices       = 0;
  mesh->num_boundary_quads = 0;
  mesh->num_quads          = 0;
  mesh->num_hexes          = 0;

  mesh->vertices      = NULL;
  mesh->quads         = NULL;
  mesh->hexes         = NULL;
  mesh->hex_adjacency = NULL;
  mesh->hex_facets    = NULL;
}

error_code mesh_load(mesh *mesh, const char *path) {
  error_code status = SUCCESS;

  mesh_init(mesh);

  FILE *file = fopen(path, "r");
  if (!file)
    return IO_ERROR;

  char buffer[256];
  while (fgets(buffer, sizeof(buffer), file) && status == SUCCESS) {
    char *it = buffer;
    while (*it && isspace(*it)) it++;

    const char vertices_cmd[] = "Vertices";
    const char quads_cmd[]    = "Quads";
    const char quads_cmd2[]   = "Quadrilaterals";
    const char hexes_cmd[]    = "Hexahedra";

    if (memcmp(it, vertices_cmd, sizeof(vertices_cmd)-1) == 0)
      status = read_vertices(mesh, file);
    else if (memcmp(it, quads_cmd, sizeof(quads_cmd)-1) == 0 ||
             memcmp(it, quads_cmd2, sizeof(quads_cmd2)-1) == 0)
      status = read_quads(mesh, file);
    else if (memcmp(it, hexes_cmd, sizeof(hexes_cmd)-1) == 0)
      status = read_hexes(mesh, file);
  }

  if (status == SUCCESS && !feof(file))
    status = IO_ERROR;

  if (status != SUCCESS)
    mesh_release(mesh);

  fclose(file);

  return status;
}

error_code save_mesh(mesh *mesh, const char *filename) {
  error_code status = SUCCESS;
  FILE *out = fopen(filename, "w");
  if (!out)
    return IO_ERROR;
  fprintf(out, "MeshVersionFormatted 2\n");
  fprintf(out, "Dimension 3\n");
  fprintf(out, "Vertices\n");
  fprintf(out, "%u\n", mesh->num_vertices);
  for (vertex_index i = 0; i < mesh->num_vertices; i++) {
    fprintf(out, "%g %g %g 0\n",
            mesh->vertices[4*i + 0], mesh->vertices[4*i + 1],
            mesh->vertices[4*i + 2]);
  }
  fprintf(out, "Hexahedra\n");
  fprintf(out, "%u\n", mesh->num_hexes);
  for (hex_index i = 0; i < mesh->num_hexes; i++) {
    for (hex_vertex_index j = 0; j < 8; j++)
      fprintf(out, "%u ", mesh->hexes[8*i + j] + 1);
    fprintf(out, "0\n");
  }
  fprintf(out, "End\n");
  fclose(out);
  return status;
}

error_code mesh_reserve(mesh *mesh, uint32_t num_hexes) {
  uint32_t num_interior_quads = 3*num_hexes - mesh->num_boundary_quads/2;
  uint32_t num_quads = mesh->num_boundary_quads + num_interior_quads;

  size_t num_quads_to_copy = num_quads;
  if (num_quads_to_copy > mesh->num_quads) num_quads_to_copy = mesh->num_quads;

  size_t num_hexes_to_copy = num_hexes;
  if (num_hexes_to_copy > mesh->num_hexes) num_hexes_to_copy = mesh->num_hexes;

  vertex_index *quads;
  if (posix_memalign((void**)&quads, ALIGNMENT,
                     4*sizeof(*quads)*num_quads) != 0) {
    return OUT_OF_MEMORY;
  }

  if (mesh->quads) {
    memcpy(quads, mesh->quads, 4 * num_quads_to_copy * sizeof(*quads));
    free(mesh->quads);
  }

  mesh->quads = quads;

  vertex_index *hexes;
  if (posix_memalign((void**)&hexes, ALIGNMENT,
                     8*sizeof(*hexes)*num_hexes) != 0) {
    return OUT_OF_MEMORY;
  }

  if (mesh->hexes) {
    memcpy(hexes, mesh->hexes, 8 * num_hexes_to_copy * sizeof(*hexes));
    free(mesh->hexes);
  }

  mesh->hexes = hexes;

  quad_index *hex_facets;
  if (posix_memalign((void**)&hex_facets, ALIGNMENT,
                     6*sizeof(*hex_facets)*num_hexes) != 0) {
    return OUT_OF_MEMORY;
  }

  if (mesh->hex_facets) {
    memcpy(hex_facets, mesh->hex_facets, 6 * num_hexes_to_copy * sizeof(*hex_facets));
    free(mesh->hex_facets);
  }

  mesh->hex_facets = hex_facets;

  hex_index *hex_adjacency;
  if (posix_memalign((void**)&hex_adjacency, ALIGNMENT,
                     6*sizeof(*hex_adjacency)*num_hexes) != 0) {
    return OUT_OF_MEMORY;
  }

  if (mesh->hex_adjacency) {
    memcpy(hex_adjacency, mesh->hex_adjacency,
           6 * num_hexes_to_copy * sizeof(*hex_adjacency));
    free(mesh->hex_adjacency);
  }

  mesh->hex_adjacency = hex_adjacency;

  mesh->num_quads = num_quads;
  mesh->num_hexes = num_hexes;

  return SUCCESS;
}

void mesh_release(mesh *mesh) {
  if (mesh->hex_facets) free(mesh->hex_facets);
  if (mesh->hex_adjacency) free(mesh->hex_adjacency);
  if (mesh->hexes) free(mesh->hexes);
  if (mesh->quads) free(mesh->quads);
  if (mesh->vertices) free(mesh->vertices);
}

static error_code read_vertices(mesh *mesh, FILE *file) {
  char buffer[256];
  if (!fgets(buffer, sizeof(buffer), file)) {
    log_error("failed to read the number of vertices in the mesh");
    return IO_ERROR;
  }

  char *end = NULL;
  mesh->num_vertices = strtoul(buffer, &end, 10);
  if (!end || (*end != '\0' && *end != '\n'&& *end != '\r')) {
    log_error("failed to parse the number of vertices in the mesh");
    return PARSE_ERROR;
  }

  if (posix_memalign((void**)&mesh->vertices, ALIGNMENT,
                     4*sizeof(*mesh->vertices)*mesh->num_vertices) != 0) {
    log_error("failed to allocate memory to store vertices");
    return OUT_OF_MEMORY;
  }

  for (uint32_t i = 0; i < mesh->num_vertices; i++) {
    if (!fgets(buffer, sizeof(buffer), file)) {
      log_error("failed to read vertex #%u", i);
      return IO_ERROR;
    }

    float data[3];
    if (sscanf(buffer, "%f %f %f 0", &data[0], &data[1], &data[2]) != 3) {
      log_error("failed to parse vertex #%u", i);
      return PARSE_ERROR;
    }

    for (int j = 0; j < 3; j++) mesh->vertices[4*i + j] = data[j];
    mesh->vertices[4*i + 3] = 0.0;
  }

  return SUCCESS;
}

static error_code read_quads(mesh *mesh, FILE *file) {
  char buffer[256];
  if (!fgets(buffer, sizeof(buffer), file))
    return IO_ERROR;

  char *end = NULL;
  mesh->num_boundary_quads = strtoul(buffer, &end, 10);
  mesh->num_quads = mesh->num_boundary_quads;
  if (!end || (*end != '\0' && *end != '\n'&& *end != '\r')) {
    log_error("failed to parse the number of quadrangles in the mesh");
    return PARSE_ERROR;
  }

  vertex_index num_vertices = 0;

  if (posix_memalign((void**)&mesh->quads, ALIGNMENT,
                     4*sizeof(*mesh->quads)*mesh->num_quads) != 0) {
    log_error("failed to allocate memory to store quadrangles");
    return OUT_OF_MEMORY;
  }

  for (uint32_t i = 0; i < mesh->num_quads; i++) {
    if (!fgets(buffer, sizeof(buffer), file)) {
      log_error("failed to read quadrangle #%u", i);
      return IO_ERROR;
    }

    uint32_t quad[4];
    if (sscanf(buffer, "%u %u %u %u 0",
               quad + 0, quad + 1, quad + 2, quad + 3) != 4) {
      log_error("failed to parse quadrangle #%u", i);
      return PARSE_ERROR;
    }

    for (quad_vertex_index j = 0; j < 4; j++) {
      if (quad[j] == 0 || quad[j] - 1 > UINT8_MAX) {
        log_error("quadrangle #%u: vertex %u is out of bounds", i, quad[j]);
        return PARSE_ERROR;
      }

      mesh->quads[4*i + j] = quad[j] - 1;
      if (num_vertices < quad[j]) num_vertices = quad[j];
    }
  }

  return reorient_quads(mesh->quads, mesh->num_quads, num_vertices);
}

static error_code read_hexes(mesh *mesh, FILE *file) {
  char buffer[256];
  if (!fgets(buffer, sizeof(buffer), file)) {
    log_error("failed to read the number of hexahedra in the mesh");
    return IO_ERROR;
  }

  char *end = NULL;
  mesh->num_hexes = strtoul(buffer, &end, 10);
  if (!end || (*end != '\0' && *end != '\n')) {
    log_error("failed to parse the number of hexahedra in the mesh");
    return PARSE_ERROR;
  }

  if (mesh->num_hexes > UINT8_MAX + 1) {
    log_error("requested number of hexahedra too high: %u", mesh->num_hexes);
    return PARSE_ERROR;
  }

  if (posix_memalign((void**)&mesh->hexes, ALIGNMENT,
                     8*sizeof(*mesh->hexes)*mesh->num_hexes) != 0) {
    log_error("failed to allocate memory to store hexahedra");
    return OUT_OF_MEMORY;
  }

  if (posix_memalign((void**)&mesh->hex_adjacency, ALIGNMENT,
                     6*sizeof(*mesh->hex_adjacency)*mesh->num_hexes) != 0) {
    log_error("failed to allocate adjacency table for hexahedra");
    return OUT_OF_MEMORY;
  }
  memset(mesh->hex_adjacency, 0,
         6*sizeof(*mesh->hex_adjacency)*mesh->num_hexes);

  if (posix_memalign((void**)&mesh->hex_facets, ALIGNMENT,
                     6*sizeof(*mesh->hex_facets)*mesh->num_hexes) != 0) {
    return OUT_OF_MEMORY;
  }
  memset(mesh->hex_facets, 0, 6*sizeof(*mesh->hex_facets)*mesh->num_hexes);

  for (uint32_t i = 0; i < mesh->num_hexes; i++) {
    if (!fgets(buffer, sizeof(buffer), file)) {
      log_error("failed to read line for hexahedron #%u", i);
      return IO_ERROR;
    }

    uint32_t hex[8];
    if (sscanf(buffer, "%u %u %u %u %u %u %u %u 0",
               hex + 0, hex + 1, hex + 2, hex + 3,
               hex + 4, hex + 5, hex + 6, hex + 7) != 8) {
      log_error("failed to parse hexahedron #%u", i);
      return PARSE_ERROR;
    }

    for (hex_vertex_index j = 0; j < 8; j++) {
      if (hex[j] == 0 || hex[j] - 1 > UINT8_MAX) {
        log_error("hexahedron #%u: vertex %u is out of bounds", i, hex[j]);
        return PARSE_ERROR;
      }

      mesh->hexes[8*i + j] = hex[j] - 1;
    }
  }

  return SUCCESS;
}

uint32_t mesh_num_boundary_vertices(const mesh *mesh) {
  uint32_t n = 0;
  for (uint32_t i = 0; i < 4*mesh->num_boundary_quads; i++) {
    if (n < mesh->quads[i]) n = mesh->quads[i];
  }

  return n + 1;
}

static bool quad_contains_edge(vertex_index *quad,
                               vertex_index a, vertex_index b) {
  for (quad_vertex_index i = 0; i < 4; i++) {
    if (quad[i] == a && quad[(i+1) % 4] == b)
      return true;
  }

  return false;
}

static bool quad_should_flip(quad_index *other_side,
                             vertex_index *oriented_quads,
                             vertex_index *quad,
                             vertex_index num_vertices) {
  for (quad_vertex_index i = 0; i < 4; i++) {
    vertex_index a = quad[i];
    vertex_index b = quad[(i+1) % 4];

    quad_index other = other_side[diagonal_index(num_vertices, a, b)];
    if (other == NO_QUAD) continue;

    if (quad_contains_edge(oriented_quads + 4 * other, a, b))
      return true;
  }

  return false;
}

static quad_index find_next_quad(quad_index *other_side,
                                 vertex_index *quads,
                                 quad_index num_quads,
                                 vertex_index num_vertices) {
  for (quad_index i = 0; i < num_quads; i++) {
    for (quad_vertex_index j = 0; j < 4; j++) {
      vertex_index a = quads[4 * i + j];
      vertex_index b = quads[4 * i + ((j+1) % 4)];

      quad_index other = other_side[diagonal_index(num_vertices, a, b)];
      if (other != NO_QUAD)
        return i;
    }
  }

  return 0;
}

static error_code reorient_quads(vertex_index *quads,
                                 quad_index num_quads,
                                 vertex_index num_vertices) {
  /**
   * Naive version without computing the adjacency, given the small number of
   * vertices that we can deal with in the first place.
   */

  quad_index *other_side;
  if (posix_memalign((void**)&other_side, ALIGNMENT,
                     sizeof(*other_side) * (size_t)num_vertices *
                     (size_t)num_vertices) != 0) {
    log_error("failed to allocate memory to compute the orientation of boundary quads");
    return OUT_OF_MEMORY;
  }

  for (size_t i = 0; i < (size_t)num_vertices * (size_t)num_vertices; i++)
    other_side[i] = NO_QUAD;

  for (quad_index num_oriented = 0; num_oriented < num_quads; num_oriented++) {
    quad_index i = find_next_quad(other_side, quads + 4 * num_oriented,
                                  num_quads - num_oriented, num_vertices);

    vertex_index *quad = quads + 4 * num_oriented;

    vertex_index tmp[4];
    memcpy(tmp, quads + 4 * (num_oriented + i), sizeof(tmp));
    memmove(quads + 4 * (num_oriented + i), quad, sizeof(tmp));
    memcpy(quad, tmp, sizeof(tmp));

    if (quad_should_flip(other_side, quads, quad, num_vertices)) {
      vertex_index tmp = quad[1];
      quad[1] = quad[3];
      quad[3] = tmp;

      if (quad_should_flip(other_side, quads, quad, num_vertices)) {
        free(other_side);
        log_error("input quadrangulation is not orientable");
        return UNORIENTABLE_ERROR;
      }
    }

    for (quad_vertex_index j = 0; j < 4; j++) {
      vertex_index a = quad[j];
      vertex_index b = quad[(j+1) % 4];
      other_side[diagonal_index(num_vertices, a, b)] = num_oriented;
    }
  }

  free(other_side);

  return SUCCESS;
}
